package topOpt;

import inf.jlinalg.*;

public class TopologyOptimization {
	//Attributes
	
	//Methods
	
//	public static void top(int nelx, int nely, double volfrac, double penal, double rmin){
//		//initialize matrix x
//		IMatrix x = new Array2DMatrix(nelx, nely);
//		IMatrix x_old = new Array2DMatrix(nelx, nely);
//		
//		//set element in x = volfrac
//		for (int i = 0; i < nelx; i++){
//			for (int j = 0; j < nely; j++){
//				x.set(i, j, volfrac);
//			}
//		}
//		
//		int loop = 0;
//		double change = 1.0;
//		
//		//Start iteration
//		while (change > 0.01){
//			loop += 1;
//			
//			//save current x to x_old
//			x.assign(x_old);
//			return;
//			
//		}
//		
//		//print
//		System.out.println(MatrixFormat.format(x));
//		System.out.println(MatrixFormat.format(x_old));
//	}
	
	public static void main(String[] args){
		//public static void top(int nelx, int nely, double volfrac, double penal, double rmin){
		int nelx = 4;
		int nely = 5;
		double volfrac = 0.8;
		double penal = 0;
		double rmin = 0;
		
		//initialize matrix x
		IMatrix x = new Array2DMatrix(nelx, nely);
		IMatrix x_old = new Array2DMatrix(nelx, nely);
		
		//set element in x = volfrac
		for (int i = 0; i < nelx; i++){
			for (int j = 0; j < nely; j++){
				x.set(i, j, volfrac);
			}
		}
		
		int loop = 0;
		double change = 1.0;
		double lk = 0;
		//Start iteration
		while (change > 0.01){
			loop += 1;
			
			//save current x to x_old
			x_old.assign(x);
			
			//FE Analysis
			//Matrix U : nelx, nely, x, penal
			//Matrix K
			//double KE = lk;
			
			//KE = lk
			double E = 1.0;
			double nu = 0.3;
			double[] k = {1/2 - nu/6, 1/8 + nu/8, -1/4 - nu/12, -1/8 + 3*nu/8, -1/4 + nu/12, -1/8 - nu/8, nu/6, 1/8 - 3*nu/8};
			IMatrix KE = new Array2DMatrix(k.length,k.length);
			double value_mul = E/(1 -  Math.pow(nu, 2));
			//set manually
			double temp[][] = {
				      { value_mul*k[0], value_mul*k[1], value_mul*k[2], value_mul*k[3], value_mul*k[4], value_mul*k[5], value_mul*k[6], value_mul*k[7] },
				      { value_mul*k[1], value_mul*k[0], value_mul*k[7], value_mul*k[6], value_mul*k[5], value_mul*k[4], value_mul*k[3], value_mul*k[2] },
				      { value_mul*k[2], value_mul*k[7], value_mul*k[0], value_mul*k[4], value_mul*k[6], value_mul*k[3], value_mul*k[4], value_mul*k[1] },
				      { value_mul*k[3], value_mul*k[6], value_mul*k[5], value_mul*k[0], value_mul*k[7], value_mul*k[2], value_mul*k[1], value_mul*k[4] },
				      { value_mul*k[4], value_mul*k[5], value_mul*k[6], value_mul*k[7], value_mul*k[0], value_mul*k[1], value_mul*k[2], value_mul*k[3] },
				      { value_mul*k[5], value_mul*k[4], value_mul*k[3], value_mul*k[2], value_mul*k[1], value_mul*k[0], value_mul*k[7], value_mul*k[6] },
				      { value_mul*k[6], value_mul*k[3], value_mul*k[4], value_mul*k[1], value_mul*k[2], value_mul*k[7], value_mul*k[0], value_mul*k[5] },
				      { value_mul*k[7], value_mul*k[2], value_mul*k[1], value_mul*k[4], value_mul*k[3], value_mul*k[6], value_mul*k[6], value_mul*k[0] }				      
				    };
			
			for (int i = 0; i < k.length; i++){
				for (int j = 0; j < k.length; j++){
					KE.set(i, j, temp[i][j]);					
				}	
			}
			
			int new_nelx = 2*(nelx+1)*(nely+1);
			int new_nely = 2*(nelx+1)*(nely+1);
			
			IMatrix K = new Array2DMatrix(new_nelx, new_nely);
			
			
			//IMatrix F = new Array2DMatrix(2*(nely+1)*(nelx+1), 1);
			double[] F = new double[new_nelx];
			for(int i = 0; i < new_nelx; i++){
				F[i] = 0.0;
			}
			
			double[] U = new double[new_nelx];
			for(int i = 0; i < new_nelx; i++){
				F[i] = 0.0;
			}
			
			//Calculate
			//double[] edof = new double[8];
			double sum = 0;
			for (int i = 0; i < nelx; i++){
				for (int j = 0; j < nely; j++){
					int n1 = (nelx+1)*(i-1) + j;
					int n2 = (nely+1)*i + j;
					int[] edof = {2*n1-1, 2*n1, 2*n2-1, 2*n2, 2*n2+1, 2*n2+2, 2*n1+1, 2*n1+2};
					for (int k = 0; k < edof.length; k++){
						sum = K.get(edof[k], edof[k]);
						//sum = sum + Math.pow(x.get(j, i), penal)*KE.get(i, j);
						
						
						K.set(edof[k], edof[k], sum);		
					}
				}
			}
			
		}
		//x_old.assign(x);
		
		//print
		System.out.println(MatrixFormat.format(x));
		System.out.println(MatrixFormat.format(x_old));
		
		//System.out.println("test");
	}
}
